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Abstract 

To describe the tunneling dynamics of a stack of two-dimensional fermionic 
superfluids in an optical potential, we derive an effective action functional 
from a path integral treatment. This effective action leads, in the saddle 
point approximation, to equations of motion for the density and the phase of 
the superfluid Fermi gas in each layer. In the strong coupling limit (where 
bosonic molecules are formed) these equations reduce to a discrete nonlinear 
Schrodinger equation, where the molecular tunneling amplitude is reduced 
for large binding energies. In the weak coupling (BCS) regime, we study 
the evolution of the stacked superfluids and derive an approximate analytical 
expression for the Josephson oscillation frequency in an external harmonic 
potential. Both in the weak and intermediate coupling regimes the detec- 
tion of the Josephson oscillations described by our path integral treatment 
constitutes experimental evidence for the fermionic superfluid regime. 



I. INTRODUCTION 



Recent experiments have demonstrated that both degenerate Fermi gases and Bose- 
Einstein condensates (BECs) can be loaded in one-dimensional optical lattices created by 
standing laser waves [1—5]. The atoms are trapped in the valleys of the periodic potential, 
and form a stack of 'pancake' shaped clouds weakly coupled to each other. When the 
laser power is large enough, the gases in the separate valleys become quasi two-dimensional 
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(quasi-2D). An additional parabolic potential, provided by external magnetic fields, and 
with corresponding oscillator length much larger than the period of the periodic potential, 
can be applied. In the ground state of this system, the fermionic and/or bosonic atoms are 
distributed in the lattice sites near the bottom of the additional parabolic trap. 

The possibility to load a BEC in a periodic potential has led to the observation of the 
Mott-Insulator phase transition [6] and the detection of Bloch oscillations and Josephson 
currents through the potential barriers separating the layers of superfluid [2-4]. To observe 
the Josephson effect, Cataliotti et al. suddenly displaced the additional parabolic potential, 
placing the stack of quasi-2D BECs out of equilibrium [2]. Josephson currents allow the 
superfluid to tunnel between different layers and perform pendulum-like oscillations around 
the equilibrium position, driven by the external harmonic potential. 

In Ref. [2], a critical Josephson current was found when the BEC was moved too far out 
of equilibrium, indicating the breakdown of superfluidity across the layers [7] . The existence 
of the Josephson currents is a direct manifestation of phase coherence across layers. For a 
one-dimensional array of dilute Fermi gases the superfluid regime is predicted to be accessible 
[8], and also in this case the Josephson effect will be a signature of superfluidity. 

In this paper, we derive an effective action, starting from the path integral representa- 
tion of the partition function of the coupled quasi-2D layers filled with two different species 
of fermions at temperature zero. From this effective action, we derive equations of mo- 
tion that allow us to study the dynamics of the phase and the density of the fermionic 
superfluid. Based on our results for the equations of motions of a fermionic superfluid in a 
one-dimensional potential, we describe what would happen if a fermionic superfluid instead 
of a BEC would be subject to the experiment of Ref. [2] as illustrated in Fig. 1. We show 
that when the interatomic interactions are weak there are loosely bound BCS-pairs that 
can tunnel coherently through the barriers that separate the potential valleys. When the 
interatomic interaction becomes stronger, confinement induced quasi-2D bosonic molecules 
with high binding energy are formed and the superfluid becomes a Bose-Einstein condensate 
of these molecules and we calculate their tunneling energy. 



In the derivation of the effective action, we follow rather closely the approach suggested 
by S. De Palo et al. in Ref. [9] and start from the path integral representation of the partition 
function for a system consisting of layers of 2D fermions 



II. THE EFFECTIVE ACTION 



Z = J Vil>l ff (x)Vip jtff (x) exp [-S y>} >(T (x), ipj,a( x )} } , 



(1) 



where the action is given by 
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Here, (3 = 1/(&bT) where T denotes the temperature and ks the Boltzmann constant. The 
three-vector notation x = (x,r) is used. The field (x) belongs to a fermion of mass m 
in layer j and spin a (| or |). The potential V ext (j) is an additional external potential that 
the fermions are subjected to. The attraction strength between the fermions is determined 
by U . The interlayer tunneling energy for a fermion is denoted by the real number ti that 
can be calculated with the approximate formula from Ref. [10] for the tunneling energy in 
a optical potential V (z) = Vq sin 2 (2-nz/X) with wave length A and depth Vq: 
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where ujl = ^8n 2 Vo/ (mA 2 ) and £l — y^l/m^L are respectively the trapping frequency and 
the oscillator length that an atom feels in in de z-direction. 

In appendix A, the reduction of (2) to an effective action is given and here we only give 
a summary. In order to grasp the most important part of the path integral in the superfluid 
state, the interaction between the fermions is decoupled by the Hubbard- Stratonovich (HS) 
transformation with the complex HS-field A^ s (x). After integration over the fermion fields, 
one is left with an effective action in terms of the HS-fields. 

However, no information about the physical density of the system can be read off from 
such an effective action. In order have access to this quantity in the effective action, we 
introduce it by multiplying the partition function with the constant 



C = J V(? s (x)V Pj (x) exp - ]T j j rf 2 x <f 5 (x) 
x Pj (x) - V>j, T (x) 4> jA (x) - (x) (x)] } . 
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Carrying out the functional integral over (j (x) alone gives S[pj (x) — ip}^ (x) ipj, 
ipj^ (x) ipj^i (x)] and thus pj (x) corresponds to the physical density of the system along any 
path. Next, the complex field Af s (x) is separated in a modulus and a phase. This phase 
is important for the low energy dynamics and therefore it is advantageous to introduce it 
explicitly 
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We then arrive at the following expression for the partition function 



Af 3 (x) VCf s {x)Ve o {x)V Po {x) exp [S eS ] , 



(5) 
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where S e s is given in appendix A, expression (A10). 



To describe the low-energy dynamics of the density and the phase of the superfluid, the 
paths along which 0j(x) and Pj(x) vary slowly in comparison to the fermionic frequencies 
(Fermi energy and binding energy) will be of importance. Along these paths, we make a 



saddle point approximation for the remaining fields 
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Af s (x) and 5(f s (x) around the saddle point values 



x) in appendix B. 
\f\x) andCfV) 
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can be treated perturbatively. The saddle point value for the effective action is calculated 
in the appendix, and given by expression (B9). The saddle point equations are 
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with n F (E) = \j{eP E + 1) the Fermi-Dirac distribution function and Ej (k) the local BCS 
energy defined by 
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Equations (8) and (9) show that the saddle point value (f*\x) can be interpreted as a 
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chemical potential Zj(x) = i(j°\x). The first saddle point equation (7) corresponds to the 
BCS gap equation, whereas the second saddle point equation leads to the BCS equation 
fixing the chemical potential Zj(x) in layer j as a function of the density Pj(x) in layer j. 

As we have introduced a momentum independent contact-interaction, equation (7) has 
to be regularized as described in [11], after which it becomes 
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T 00 (E) 
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2n F [Ej(h)] - 1 
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where Too (E) is the low-momentum limit of the T-matrix. This equation has no ultraviolet 
divergences. In two dimensions, at low energy, T o (E) is given by [12,13] 



1 



Too (E) 



m 
T 



7T 



In [E/E b ] + i 



(11) 



where E b is the energy of the 2D bound state that always exists in two dimensions. For the 
optical potential V sin 2 (2nz/X) described earlier, the binding energy of the quasi-2D bound 
state is given by [16] 
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E b = exp 
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(12) 



with a the scattering length of the fermionic atoms and C 0.915. 



III. JOSEPHSON CURRENT AT T = 



A. Equations of motion for density and phase 

We now proceed with an analysis at T = and with the assumption that the energy ti 
is small compared to the other energies, such that a perturbational expansion with ti as a 
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small parameter is possible. In this case, equations (8) and (10) can be solved analytically 



for 



Af\x) and Zj(x) (see also [12]) : 

Zj(x) 



l27T Pj (x) 



V m 

n Pj {x) 



m 
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As we want to study the current perpendicular to the layers in which the atoms are 
confined, we have to calculate the terms in the effective action that couple the different layers. 
In our perturbational expansion of the effective action (B9), the lowest order contribution 
comes from the term in the self energy (B7) proportional to t\. We find that the contribution 
equals 



3 ^ w„=(2n+l)7r//3 



d 2 k 



(27T) 



E] +1 (k) ul + E](k) 
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cos [0 j+ i(x) -6j(x)]}, 



(15) 
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with £j (k) = k 2 1 (2m) — Zj(x) the free fermion dispersion and Ej the BCS dispersion relation 
(9). In this expression, the part proportional to cos[0j + i(x) — Oj(x)] is responsible for the 
Josephson tunneling between the layers and we write it symbolically as 



Stunnel = ~Y1 J Q d T J dx T j+1J COS [9j+l(x) ~ 6j{x)\ 



(17) 



Evaluating this term with the supposition that the gap 
Zj (x) vary slowly with j leads to 
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A similar contribution to the energy can be obtained in a BCS-approach, following Ref. [14]. 
We now take as an approximation for the effective action 



Sj [Bj{x),Pj{x)\ = Sl% + Stunnel, 



(19) 



with Slff given by (B9) from which 



A°(x) 



and Cj{ x ) have been eliminated using equations 



(13), (14). Having obtained an expression for the action that only depends on 9j(x) and Pj(x) 
we can finally proceed with deriving equations of motion by extremizing Sj with respect to 
these fields. Because we want to give a dynamical interpretation to these equations, we write 
them in real time (id T — > d t ). Extremizing with respect to the phase field 6j(x) results in 



Pj ( x ) V0j (x) ■ Vpj (x) 



Am 



+Tj,3-ism [6j(x) - Oj-^x)] 
-T j+1J am [9 j+ i(x) - 8j(x)] , 



(20) 
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and the derivative with respect to pj (x) yields 



-dt- 1 ^- = - — 3 + V ext (?) +Zi-u 



cos [9j + i(x) — Oj(x)] 



dpj (x) 

cos^x) - O^x)} . (21) 



We have calculated the derivative dTj+ij/dpj if the density varies smoothly with the layer 
index (pj+i ~ Pj ~ Pj-i) and is constant in the plane: 



dpj(x) 2 [27rp j (x)/m + E b ] 
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(22) 



B. Oscillations of the superfluid in an optical lattice 

We now introduce the wave function 



m = f-4^^ (t \ (23) 

with pj and 9j only depending on time and layer index j. That is, we assume that within 
a layer, the density and the phase are homogeneous, but that they can still vary over time 
and over layers. The wave function (23) obeys the Schrodinger equation 



-^' w F lT /m / * i cos - + cos ^ - • (24) 

+ Zirpj/m) 

where we can take in the approximation of slowly varying density that we used before, namely 
ipj-i ~ ipje 1 ^- 1-6 ^ and ipj + i ps ipje % ^ j+1 ~ 6j \ so that we have in the limit rcpj/m <C E b that 

*^ = (2^ (j) + 2 Zj - 2p) ^ - 2 ^ + Eb (^-i + ^+i) , (25) 

which is the nonlinear discrete Schrodinger equation for an array of Bose-Einstein conden- 
sates [15], where the boson tunneling matrix element is given by t\j {2npj/m + Eb). This is 
a decreasing function of Eb, which has the physical consequence that for increasing molecular 
binding energy the coupling between the layers becomes weaker and that in the Bose-Einstein 
limit {E h — > oo), no tunneling of molecules occurs. This may be related to the fact that the 
binding together of atoms in molecules is strongly influenced by the the confinement poten- 
tial and the quasi-2D nature of the gas (as can be seen from eq. 12). In the intermediate 
states in the calculation of the amplitude of the tunneling process, the molecular state is 



6 



unlikely to survive as a bound state. Hence, the molecular binding energy is to be added to 
the energy barrier for tunneling. 

We know from [2], that in the bosonic limit equation (25) allows for oscillations where 
the phase difference 9 j+1 — 9j is locked to a constant value A6> and the differential equation 
governing the dynamics of the center of mass coordinate and this phase difference A9 is of 
the pendulum type. Our numerical simulations of (24) suggest that in the BCS-regime there 
is a similar current through the array. 

iFrom equations (20) and (21), we can derive a simplified analytical formula to estimate 
the oscillation frequency in an external harmonic potential V ext (j) = Qj 2 . Assuming that 
the phase difference between neighboring layers is constant, this becomes 

d t R = 2 (Tjj-i) sin (9j - 9^) , (26) 

where 

is the average of over the lattice sites. The initial density profile pj(t = 0) is calculated 
from the chemical potentials Zj through (14) and these chemical potentials are derived from 
(21), Zj — ii — V ext (j), for a particular choice of /i. The last two terms from equation (24) 
seem to have only a minor influence on the frequency so that we omit them in the analytic 
calculation. Equation (21) then becomes 

d t (9j - Oj-!) = -AVLR. (27) 

For small oscillations, (26) and (27) lead to the frequency 



CO 



= Jmr~T). (28) 



Taking 40 K atoms with a central density of n = 10 9 cm~ 2 , an optical wavelength A = 754 
nm, and an axial frequency u a = 2n x 24 Hz, we plot in Fig. 2 the analytical frequency 
(curves) obtained from equation (28) and compare it with a numerical calculation (symbols) 
based on equation (24). The oscillation frequency is plotted as a function of the inverse 
scattering length which appears in expression (12) of the binding energy. Fig. 2 shows that 
the estimation (28) agrees almost perfectly with the numerical results in the bosonic limit, 
and that also in the BCS regime there is reasonable quantitative agreement. The inset of Fig. 
2 shows that far in the bosonic regime (a > 0), the oscillation frequency decreases rapidly 
as an exponential function of 1/a. Nevertheless in the cross-over regime, the oscillation 
frequencies are high enough for the observation of Josephson currents to be useful as a tool 
to investigate the superfluidity of an ultracold Fermi gas, analogous to the experiments of 
Ref. [2] in the bosonic case. 



IV. CONCLUSIONS 

We have derived an effective action and the resulting equations of motion to describe 
the dynamics of a fermionic superfluid in a layered system and applied this formalism to 
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study the center of mass motion of an atomic Fermi gas in the potential formed by an op- 
tical standing wave. We find that a Fermi gas in the BCS regime can perform superfluid 
oscillations through the optical lattice, similar to those that have been observed for conden- 
sates of bosonic atoms [2-4], when the gas is not in equilibrium in the harmonic trapping 
potential superimposed on the optical lattice. An analytical approximate expression (28) 
for the oscillation frequency is derived and the predictions of this expression are tested with 
numerical simulations of the full equations of motion (24). For a superfluid Fermi gas in 
the BEC regime, we find that the tunneling is suppressed when the molecular binding en- 
ergy becomes large. We conclude that superfluidity in Fermi gases can be revealed through 
Josephson currents in optical lattices if the Fermi gas is either in the BCS regime, or in the 
weakly-bound molecular BEC regime. 
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APPENDIX A: DERIVATION OF THE EFFECTIVE ACTION 

In order to grasp the most important part of the path integral in the paired-fermion state, 
the interaction between the fermions is decoupled by the Hubbard- Stratonovich transforma- 
tion, which transforms (2) into 



lands. 




(Al) 



with 
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Af s (x) 



U 
V 2 



+ E ^l* ( x A dr ~^ m ~ + Vext 0') - A* J ^> ( x ) 



<T=±1 



- Af s (x) (*) (*) - Af s 't (*) (*) (*) 



(A2) 



In order to keep track of the total density, we introduce the constant C (see expression 
(4)). In order to investigate the BCS gap and the phase we separate the complex field 
A^ s (x) in a modulus and a phase (expression (5)) and also transform the fermion fields as 
4>j,a i x ) — ¥ ^j,a ( x ) e' 0J '^/ 2 . Additionally, we shift the field i(f s (x) according to 



iCf S (x)^i("*(x)+id. 



HS 



9 s (x) , (V9 3 (x)Y 



+ 



8m 



+ V ext (j) - n 



(A3) 



and use the Nambu spinor notation m (x) = ypj^ (x) , (x)J . After this procedure, the 
partition function can be written as 

Z oc J |Af s (x)\ Vr${x)Vr}j{x)V \dJ? s (x)\V0j{x)V(;? s {x)Vpj{x) exp {-S {2) } , (A4) 

with the action = 5*0 + [r)j~(x), r)j(x)] where 
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x < 



U 



+ 



2 8m 
does not contain the fermion fields any more and 



+ V ext 



Pi (x) 



(A5) 



5(3) = E [ dT f<?Ari]{x) 



' .V0,-(a;)„ .V 2 9Ax)' 
d r - i j> — V -i- 3 y ' 



2m 



Am 



+ (-^-<f 5 Wk3- Af 5 (x) 

+ 



<7l 



, 2m ,J ' 7 - 3 
t lV } (x) e^-^^awj+i (x) + h.c] } 



(A6) 



is the part of the action that still depends on them. The Pauli matrices are denoted by 
(7j. Since is quadratic in the fermion fields, the path integral over these fields can be 
performed, resulting in 



J P77t(x)P7fc(x)exp{-S( 3 >} = deti-G- 1 ], 



(A7) 
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where the Green's function G is a matrix in coordinate space as in layer space and is given 
by 



-G 1 (x,j; x',j') =5{x- x') ^6jj, 

/ V 2 



d T - i ^ V - i — (T 



2m 



Am 



<7l 



The partition sum can then be written as 

Z oc J V \Af s (x)\v9 j (x)VC J HS (x)Vp j (x) exp {-S cS } 

where 

s cS = s + Tt[h-g- 1 ) 

with So given by (A5) and the Green's function given by (A8). 



(A8) 



(A9) 



(A10) 



APPENDIX B: SADDLE POINT EXPANSION FOR THE EFFECTIVE ACTION 

iFrom the path integrations in (A9), we ultimately want to extract equations of motion 
for 9j(x), Pj(x). To describe the low-energy dynamics, the contributions with 9j(x), Pj(x) 
varying slow in comparison with the fermionic frequencies (Fermi energy and binding energy) 
will be important. Along the paths of slowly varying 0j(x), Pj(x), we can make the saddle 
point expansion in the fields A^ 5 , (f s , setting 



Af s (x) 



Af (x) 



+ 5 



Af s (x) 



where also 



(Bl) 
(B2) 



A^ (x) and (x) will vary slowly in comparison to the fermion frequencies. 



— izj leads to 

G- 1 = Go 1 + S, 
where the saddle point contribution is given by 

-Go 1 x',f) =8(x- x') 8 j3 - 



Expanding the Green's function (A8) around the saddle point values A^ and (, 



.(o) 



V_2 

2m 



d T a + [ - — - Zj o- 3 - 



A {0) 



(B3) 



(B4) 



This saddle point contribution can be diagonalized by going to Fourier space : 



Gq 1 (k,w, j;k',u/, j') = Sjj'Su^S (k- k') 



iua -{*- *,) a 3 + 



0"! 



(B5) 
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and thus 



Tr 



In 



-Go 1 



x In 



w„=(2n+l)7r/ / 9 ' 



(2vr) 2 



k 2 /(2m) - ^-(x) 



In this expression, the x dependence of Zj(x) and Af\x) has been reintroduced, still 
assuming that these fields vary slowly in comparison to the relevant fermion frequencies. 
The self energy E equals 

-E (xj; x',f) =5(x- x') S jf \-i<j 3 5(? s (x) - (5 A" s (x) )a, 



+ 



(B6) 



- « 



V0,- .V 2 # ? 



so that 



det 



2m 4m 
+5 (x - x') [^ +1 ^ 1 e^+ l( * ) -^ )W V 3 
+<5 J _ iy t 1 e- J ^+ l(:E) -^ (a;)]<T3/2 C T3] , 

-G- 1 ] = exp {Tr [in (-G^ 1 )] } = exp {Tr [in (-G^ 1 - s)] } 



(B7) 



exp Tr [in (-Go 1 )] - £ 



Tr [(-G ^y 



n=l 



n 



(B8) 



where Tr denotes the trace over all the variables (coordinates, imaginary time, layer index 
and Nambu space). Expanding the action, S e g, expression (A10), around the saddle point 
values 



Aj > \x) and Zj(x) = i(j V) (x), and using the result (B8) to lowest order, we find 



•(0) 



^ = n[ln(-G^]+?jf dT / d2 



Af (x) 2 



U 



+ 



Pj (x) 



(B9) 



2 8m 

The saddle point values Af\x) and ^(x) = i(j°\x) are now equal to the extremum of 
(B9) and they fulfill the equations (7) and (8). The lowest order Green's function G is 
given by (B6). 
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FIGURES 



FIG. 1. Illustration of the experiment where the center of the harmonic+optical trap, confining 
quasi-2D clouds of atoms, is suddenly moved, which causes the oscillation of the superfluid. The 
grey wave represents the magnetic+optical trapping potential and the ellipsoids are the quasi-2D 
atomic clouds. 

FIG. 2. The center of mass oscillation frequency of a zero temperature fermionic superfluid in a 
one dimensional optical potential is shown as a function of the inverse scattering length for different 
values of the optical potential measured in units of the recoil energy E R . The lines are calculated 
with formula (28) for V /E R = 6 (solid line), V /E R = 8 (dashed line), V /E R = 10 (dashed 
line). The left of the curve Vq/Er = 6 is shown dotted because there the condition t\ <C |A| is not 
fulfilled anymore. The symbols show the frequency of a sinusoidal fit to the full numerical solution 
of the equations (24) for Vq/E r = 6 (squares), V$/E R = 8 (circles), Vq/Er = 10 (triangles). 
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